Introduction
With the coming of
the post-genome research era, more attention has been paid to the strategy of
model biology research. The model plants A.
thaliana and rice have been widely used in various research fields of
botany. However, A. thaliana, which
is a dicotyledon, has no biological characteristics
and related agronomic traits of cereal crops, such as tillering
and panicle differentiation, and is far away from monocotyledons (Blümel et al. 2015).
Rice (Oryza sativa), whose genome sequencing has
been completed early, is a model plant in Gramineae,
but it has evolved very early from the cold-season Gramineae,
and is far away from temperate crops.
Consequently, a novel model plant is needed to study some of the
problems associated with gramineous plants. B. distachyon
is a subfamily of Poaceae, which is closely related
to many important cereal crops such as wheat, barley, and rice (Brutnell et al.
2015). In addition, compared to rice, whose genome is around 430 Mb (Sasaki and Antonio 2004),
the genome of B. distachyon is just
approximate 272 Mb (Initiative 2010). B. distachyon, as a model plant
for grass crops, according to its biological characteristics of short life
cycle, and easy to cultivate and genetically transform, could offer a suitable
and meaningful foundation for many aspects of plant biology research (Scholthof et al.
2018).
Drought is one of the abiotic stresses that seriously affect plant
growth and development, and is also the main limiting factor for global crop
distribution and yield increase (Farooq et al.
2009; Matiu et al. 2017).
Plants have evolved sophisticated mechanisms to respond to and adapt to drought
stress at the molecular and physiological levels. However, previous molecular
studies about stress response mainly concentrated on genes with protein-coding
function (Sakuma et al. 2006; Huang et al. 2008).
Recently, non-protein-encoding transcripts have also proven to be
important regulators of gene expression (An et al. 2018; Nejat and Mantri 2018).
Long non-coding RNA (lncRNA) is a kind of transcript
with a length of more than 200 nucleotide (nt),
lacking coding sequence (CDS) or open reading frame (ORF), and no or almost no
protein coding potential (Do et al. 2019). Increasing evidences show that plant lncRNAs are critical regulators for growth and development
as well as response to stresses. LncRNAs could
regulate the flowering time (Heo and Sung 2011),
pollen development (Huang et al. 2018), fruit senescence (Wang et al. 2018), photomorphogenesis (Wang et al. 2014b), as well as biotic
and abiotic stress responses (Nejat and Mantri
2018). Functional analyses of lncRNAs suggest
that they take part in gene expression regulation not only at transcriptional
level, but also at post-transcriptional level through many complicated
mechanisms, such as precursors for microRNAs (miRNAs)
and other small RNAs or as target mimetics of miRNAs (Qin et al. 2017). At present, most
studies on drought stress response in the main gramineous
crops still focused on physiological levels or protein-coding functional genes,
however, little is known about how non-coding transcripts exert an influence on
the response and regulation of drought stress.
In this research, B. distachyon, the optimal model organism for gramineous crops, was taken as experimental material, and
the high-throughput sequencing technology was employed to sequence the transcriptome of drought stress and normal growth plants,
respectively. The drought-responsive lncRNAs were
identified through the comparative transcriptome
analyses to provide a theoretical foundation and new insight for understanding
and improving the molecular regulation mechanism of drought response of gramineous crops.
Materials and Methods
Plant materials, growth conditions and drought stress treatment
Cultivar of B. distachyon Bd-21 was used as the
experimental material. After removing the inferior palea,
seeds were placed in a clean petri dish covered with moist filter paper to
germinate. One day later, the petri dish was transferred to a refrigerator at 4°C
for 8 days. Then the seeds were planted in nutrient pots in a light incubator,
with growth conditions: 25°C, light for 16 h / dark for 8 h. When the seedlings
grew to the 7th leaf stage, the drought treatment was carried out by
means of water control. The control was watered once every other day, while the
treated group grew seven days without watering. Then the leaves of the
seedlings were sampled and frozen by the liquid nitrogen, storing in an
ultra-low temperature refrigerator at -80°C for following testing.
RNA isolation and qualification
The methods of total RNA isolation and reverse
transcription refer to the description of Ma et al. (2019). The degradation
and contamination of total RNA were detected by 1% agarose
gels. The purity, concentration and integrity of RNA were checked by the Nano
Photometer® spectrophotometer (IMPLEN, CA, USA) and Qubit®
RNA Assay Kit in Qubit® 2.0 Fluorometer
(Life Technologies, CA, USA), and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA),
respectively.
Library construction,
sequencing and clustering
The construction
of cDNA library, sequencing and clustering were
conducted by the Novogene Bioinformatics
Technology Co., Ltd. in Beijing, following the procedures of relevant kits,
systems and machines.
Data processing and bioinformatics
analysis
The data from Illumina
sequencing were processed further and analyzed by the Novogene
Bioinformatics Technology Co., Ltd.
(1) Quality control (QC) was checked by raw
data, clean data, Q20, Q30 and GC content of the clean data. Raw data (raw
reads) of fastq format were firstly processed through
in-house perl scripts. In this step, clean data
(clean reads) were obtained by removing reads containing adapter, reads on
containing ploy-N and low quality reads from raw data. At the same time, Q20,
Q30 and GC content of the clean data were calculated. All the downstream
analyses based on the clean data with high quality. The correlation of gene
expression levels between samples is based on the square of the Pearson
correlation coefficient (R2) > 0.92 (ideal sampling and test conditions),
and R2 should be at least greater than 0.8 in a specific test operation.
(2)
Reference genome and gene model annotation files were
downloaded from genome website (http://www.ncbi.nlm.nih.gov/nuccore/288572679)
directly. Index of the reference genome was built using Bowtie v2.0.6 and
paired-end clean reads were aligned to the reference genome using TopHat v2.0.9.
(3)
The mapped reads of each sample were assembled by
Cufflinks (v2.1.1) (Fig. S1) (Trapnell et al. 2010). According to the
assembly results of cufflinks, the structural features of lncRNA
as well as the functional qualities of non-encoded proteins, a series of severe
screening criteria were established. Then these lncRNAs
could be selected and taken as the final candidate lncRNA
set in terms of the following five steps (Fig. S2).
(4)
Coding potential analyses included Coding
Potential Calculator (CPC) as well as Pfam-scan. The
clarification of both coding and non-coding transcripts is mainly based on CPC
(0.9-r2) to assess the extent and quality of ORFs in transcripts and search
known protein sequence databases (Kong et al. 2007). Each transcript
was translated in all three possible frames and used Pfam
Scan (v1.3) to identify occurrence of any of the known protein family domains
documented in the Pfam database (Punta et al. 2012).
If a transcript has a Pfam hit, it would be removed
in the following steps.
(5)
The function of lncRNA could be
mainly achieved by cis
or trans
acting on the target gene with protein-encoding function (trans analysis is performed when the number of samples is more than
5). The basic principle of cis acting target gene prediction is that the function of lncRNA is related to the protein-coding gene adjacent to
its coordinates, so the adjacent (upstream and downstream 10 k/100k)
protein-coding genes are screened as their target genes. The main functions of lncRNA were predicted by functional enrichment analysis of
target genes.
(6)
Quantification of gene expression level was carried out
via the Cuffdiff (v2.1.1), which was used to
calculate FPKMs of both lncRNAs and coding genes in
each sample (Trapnell et al. 2010). FPKM means fragments per kilo-base of exon
per million fragments mapped, calculated based on the length of the fragments
and reads count mapped to this fragment. The threshold of differential
expression varies from sample to sample. The transcript or gene with a P-adjust
< 0.05 is the differential expression threshold for biological replicates,
while the P-adjust <0.05 and |log2 (fold change)| > 1 is a significant
difference threshold for non-biological replicates.
(7)
Gene Ontology (GO) enrichment analysis of genes with
differential expression or lncRNA target genes were
performed through the GO seq R package, in which gene
length bias could be corrected. The corrected P value of GO term which was less
than 0.05 would be considered significantly enriched by differential expressed
genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource to
recognize functions of the biological system (http://www.genome.jp/kegg/). The
KOBAS software was used to detect the statistical enrichment of differentially
expressed genes or target genes of lncRNA in KEGG
pathways.
(8)
Alternative splicing (AS) events were divided into 12
basic forms via the software Asprofile v1.0 (Florea et al. 2013).
The amount of AS events in each sample was estimated, respectively.
qRT-PCR validation
6 lncRNAs and 6 mRNAs were
chosen to detect the expression level via
qRT-PCR. The reaction system and
procedures according to the description of Feng et al. (2015). The fold
variation was calculated by 2-△△Ct method (Livak and Schmittgen 2001). The β-actin was taken as internal control.
Results
Evaluation of sequencing
quality
The Illumina HiSeqTM2500 was performed in order to obtain the
drought responsive lncRNAs in B. distachyon. In this experiment,
128,474,682 and 116,766,272 raw reads were acquired in CK and drought
treatment, respectively. The original sequencing sequences produced by
sequencing contain low-quality reads and some adapter. The raw reads must be
filtered in order to attain clean reads to assure the quality of bioinformatic analysis, because clean reads are the basis
of subsequent analysis. The clean reads acquired after filtering were
122,429,492 and 111,457,314, separately (Table S1).
The correlation of gene expression levels between samples is an
essential index to verify the reliability of the experiment and the rationality
of sample selection. The closer the correlation coefficient is to 1, the higher
the similarity in expression patterns between samples (Fig. S3).
Assembly and screening of lncRNA
Cufflinks (Trapnell et al.
2010) used the results of the TopHat
alignment to assemble transcripts, to obtain as small a collection of
transcripts as possible, and to quantify transcripts. Based on the splicing
results of cufflinks, the candidate lncRNAs were
determined by a 5-step screening method. Fig. 1 shows the statistics of the
number of transcripts produced by step1-step5 during the basic screening
process of lncRNAs.
The candidate lncRNAs were analyzed by CPC
and Pfam-scan, and the
non-coding transcripts identified by the above analyses were counted. Fig. 2
visually displays the common and unique numbers of the predicted non-coding
transcripts. The number of non-coding transcripts identified by CPC and Pfam-scan was 28,956 and 29,550, respectively. The common
part of these prediction results (28,436) was taken as the lncRNAs
data set predicted by this analysis for subsequent analysis.
Alignment analyses of
expression level
By comparing the
box plot (Fig. 3A) and density profile (Fig. 3B) of the quantitative results
between CK and drought, the differences in the overall expression level of the
samples under different treatments were visually reflected. Currently, among
the RNA-seq technologies, FPKM is the most popular
measure for evaluating gene expression levels (Trapnell et al. 2010). The FPKM values
obtained by cuffdiff quantification could be used to
compare the expression levels of mRNA and lncRNA. The
results showed no molecular type preference (Table S2). The FPKM values of lncRNA and mRNA are shown in Table S3.
The expression values of lncRNA and mRNA transcripts
of all samples were averaged, and the violin plot (Fig. 4) was drawn with the
log10 (RPKM+1) values to visually reflect the difference in the whole
expression level of lncRNA and mRNA. The results
indicated that the expression level of lncRNA was
lower than that of mRNA. Differential transcript analysis was performed using cuffdiff.
Using q-value<0.05 as the threshold of
screening, the volcano plot can be used to visually see the whole distribution
of differential transcripts or genes (Fig. 5). The red dot and blue dot in the
figure indicate the up-regulated (2,497) and down-regulated (973) transcripts
with significant differential expression, respectively (Fig. 5).
Functional prediction of lncRNA target gene
The (upstream and
downstream 100 k) protein-coding genes adjoined lncRNAs
were screened as their target genes, and the main functions of lncRNAs were predicted by functional enrichment analysis of
target genes. The predicted results of cis target genes are shown in Table S4. The GO classification
plot (Fig. 6A and Fig. 7A) represented differentially expresses lncRNA target genes and mRNAs, which intuitively reflected
the number distribution of differential genes in GO term enriched in biological
processes, cellular components, and molecular functions.
Fig. 1: Screening statistics of lncRNA.
The x-axis is the screening step, and the y-axis is the number of transcripts
after the corresponding step is filtered
Fig. 2: Venn diagram shows number of lncRNAs
of predicted non-coding transcripts
Fig. 3: Comparison of gene expression levels under different
treatments. A, FPKM box plot, the
x-axis is the sample name, the y-axis is log10 (FPKM+1), and the box plot for
each region corresponds to five statistics (top to bottom, maximum, upper
quartile, median, lower quartile and minimum); B, FPKM density distribution profile, the x-axis is log10 (FPKM+1),
and the y-axis is the density of the gene
In vivo, distinct genes collaborate with each other to exert their biological
functions. Pathway significant enrichment analysis e mploys
a hypergeometric test to ensure pathway which was
significantly enriched in specific genes compared to the whole genome
background. The scatter plots represent the results of the KEGG enrichment
analyses (Fig. 6B and Fig. 7B). The degree of KEGG enrichment is measured by
the rich factor, the q-value, and the number of genes enriched in this pathway.
Rich factor refers to the ratio of the number of genes in the pathway entry in
the lncRNA target genes with differential expression
to the total number of genes in the pathway. The larger the rich factor, the
higher the degree of enrichment. Q-value is the p-value after multiple
hypothesis test correction. The value of q-value is [0, 1], the closer to zero,
the more significant the enrichment. We selected the 20 most significant
pathway entries for enrichment here (Fig. 6B and Fig. 7B). According to the
results of KEGG enrichment analyses, clustering analysis was performed on the
significantly enriched KEGG pathway. In the analysis process, we selected the
most significant 20 metabolic pathways to analyze the difference in expression
(Fig. 6C and Fig. 7C). The expression of lncRNA in
the photosynthesis-antenna proteins pathway in CK was higher than that in
drought treatment, while the expression of lncRNA in
Glycolysis/Gluconeogenesis, Pentose and glucuronate interconversions, Glycerolipid
metabolism pathways was lower than that in drought treatment. The mRNA
expression levels in CK were higher than drought treatment in
photosynthesis-antenna proteins, metabolic pathways, carbon fixation in
photosynthetic organisms, and lower than drought treatment in pyruvate
metabolism, ascorbate and aldarate
metabolism, degradation of aromatic compounds pathways.
The response of plants to drought stress is a complex
regulatory system including multiple genes and
Fig. 4: Comparison of lncRNA and mRNA
expression levels under drought treatment. A,
violin plot of lncRNA and mRNA expression, the x-axis
is the gene name, the y-axis is log10 (FPKM+1), the width of each violin
represents the number of points at that expression level; B, FPKM density distribution profile of lncRNA
and mRNA, the x-axis is log10 (FPKM+1), and the y-axis is the density of the
gene
Fig. 5: Volcano plot of
differential expression transcript under drought treatment.
Significantly differentially expressed transcripts are represented by red dots
(upregulated expression transcripts) and green dots (downregulated expression of transcripts); x-axis represents
fold change in transcripts in different samples; y-axis represents statistical
significance of differences in transcript expression changes
regulatory networks. According to the functional annotations of predicted target
genes, we divided them into five categories, including three
photosynthesis-related genes, two transcription factors, one peptide, five
ribosomal proteins, and one glycoside hydrolase (Table S5). A total of 90 lncRNAs with differential expression were involved. Except
that the transcription factor wox11 was up-regulated under drought treatment,
the remaining target genes were down-regulated (Table S5).
The regulation mechanism of lncRNA on its
target gene is very complicated. The results of this study showed that lncRNA regulated its target genes mainly in two ways,
one-to-many and many-to-one. The former refers to that one lncRNA
could regulate multiple target genes (Table S6), while the latter, in contrast,
means that one target gene could be regulated by multiple lncRNAs
(Table S7).
AS analysis
The same gene may
have different AS patterns, which greatly increase the ability to encode genes.
Different AS patterns allow the same gene to produce multiple different mature
mRNAs, ultimately producing different proteins. AS is an important pattern by which
lncRNA regulates gene expression. The lncRNA could form a complementary double strand with the
transcript of protein-encoding gene, which interferes with the splicing of mRNA
and then could produce different AS events (Bardou et al. 2014; Gonzalez et al. 2015). The gene model
predicted by Cufflinks (Trapnell et al. 2010) was employed to
classify and quantify the AS events of each sample using AS profile (Florea et al. 2013)
software. 12 categories of AS events include (1) TSS: Alternative 5’ first exon
(transcription start site), (2) TTS: Alternative 3’ last exon (transcription
terminal site), (3) SKIP: Skipped exon (SKIP_ON,SKIP_OFF pair), (4) XSKIP:
Approximate SKIP (XSKIP_ON,XSKIP_OFF pair), (5) MSKIP: Multi-exon SKIP
(MSKIP_ON,MSKIP_OFF pair), (6) XMSKIP: Approximate MSKIP (XMSKIP_ON,XMSKIP_OFF
pair), (7) IR: Intron retention (IR_ON, IR_OFF pair), (8) XIR: Approximate IR
(XIR_ON, XIR_OFF pair), (9) MIR: Multi-IR (MIR_ON, MIR_OFF pair), (10) XMIR:
Approximate MIR (XMIR_ON, XMIR_OFF pair), (11) AE: Alternative exon ends (5’,
3’, or both), (12) XAE: Approximate AE. The number distribution of AS is
similar in the two treatments, with the largest number of TSS and TTS and the
least number of XMSKIPs (Fig. S4).
Fig. 6: Functional enrichment of differentially expressed lncRNA target gene under drought treatment. A, GO classification plot; B, scatter plot of the KEGG enrichment
analysis. The top 20 enriched target genes of differentially expressed lncRNAs. The y-axis represents the
pathway name and the x-axis represents the Rich factor. The size of the point
indicates the number of differentially expressed genes in the pathway. The
color of the point corresponds to a different q-value range; C, clustering analysis of the significantly
enriched KEGG pathway. Highly expressed metabolic pathways were indicated in
red, and relatively low metabolic pathways were indicated in green. Within each
parenthesis after the metabolic pathway tag is the number of genes
Validation of lncRNAs and mRNA
expression
6 lncRNAs and 6 mRNAs, with differential expression, were
chosen randomly in order to validate the reliability of high-throughput
sequencing by qRT-PCR.
These results showed the similar expression patterns contrasted with the transcriptome profile (Fig. 8, Table S2), indicating that
the transcriptome sequencing results were true and
reliable.
Discussion
Drought is one of the severe abiotic stresses that
influence plant growth and metabolism in various ways. For crops, the damage to
photosynthetic systems caused by drought stress is the main cause of yield
reduction (Bodner et al. 2015). Photosynthesis includes photoreaction and
dark reaction. The photosynthetic electron transfer chain in photoreaction is
mainly composed of four complexes on the thylakoid membrane: Photosynthetic
System II (PSII), Cytochrome b6f complex (Cytb6f), Photosynthetic System I (PSI)
and ATP synthase complex (Yadav et al. 2017). The Cytb6f
mediates the electron transport from PSII to PSI, which is the primary part of
the proton gradient formed on both sides of the thylakoid membrane. It is
critical in regulating the electron flow direction of the electron transport
chain as well as the state transition between PSII and PSI. In addition, the
Cytb6f is also the regulatory center of photosynthetic phosphorylation (Finazzi et al. 2016).
Ferredoxin (Fd)
is a soluble protein that could accept electrons from PSI and contribute
electrons to Fd: NADPH reductase
(FNR), which produces NADPH for the Calvin cycle and ultimately generates ATP
as well as releases oxygen (Zanetti and Allverti
2018). Two doubled haploid lines of winter triticale were treated with
water stress by Hura et al. (2019). The results showed that in both lines, water
stress caused a decrease in the content of Cytb6f protein, and the quantum
efficiency of electron transport from the original electron acceptor to the
final receptor was correspondingly reduced. In photosynthesis, excessive
reactive oxygen species (ROS), produced by reduced protein levels of the
Cytb6f, would damage the components and structures of the photosynthetic
system, and ultimately limit plant growth and yield (Yamori et al. 2016; Hura et al. 2018). Under experimental
conditions, the expression level of the Cytb6f was decreased, possibly to
protect the PS I electron acceptor from excessive reduction. The expression
levels of the above components in photosynthetic system
Fig.7: Functional enrichment of differentially expressed mRNA
under drought treatment. A, GO
classification plot; B, scatter plot
of the KEGG enrichment analysis. The top 20 enriched mRNAs with differential
expression. The y-axis represents the pathway name and the x-axis represents
the Rich factor. The size of the point indicates the number of differentially
expressed genes in the pathway. The color of the point corresponds to a
different q-value range; C,
clustering analysis of the significantly enriched KEGG pathway. Highly expressed
metabolic pathways were indicated in red, and relatively low metabolic pathways
were indicated in green. Within each parenthesis after the metabolic pathway
tag is the number of genes
were down-regulated after drought stress treatment, which is
in line with the results obtained in this study. YP_002000536.1,
YP_002000516.1 and YP_002000499.1, which are target genes (the functions are annotated as
FD, Cytb6f, and Cyt f
large domain respectively,
regulated by the corresponding lncRNA, see Table S6),
were down-regulated after drought treatment, possibly related to the role of FD
and Cytb6f in photosynthetic electron transport. The decreased expression under
stress can reduce electron transfer to limit the production of ROS and avoid
structural damage to protect the photosynthetic system.
Fig. 8: Relative expression of 12 selected genes measured by qRT-PCR under drought treatment. 6 mRNAs: BRADI1G47820,
BRADI5G25810, BRADI1G46602, BRADI4G37860, BRADI3G39467, BRADI2G11730; 6 lncRNAs: XLOC_253993, XLOC_168454, XLOC_301959,
XLOC_308502, XLOC_253993, XLOC_304290
A large part of the plant genome
is used for transcriptional regulation, for instance, 57 MADS-box genes have
been identified in B. distachyon,
involving in abiotic stress responses, including salt, drought, and
low-temperature (Wei et al. 2014). The expression patterns of these TFs were
different, among of which, 15 were up-regulated while 2 were down-regulated.
APETALA1/FRUITFULL belongs to a family of the MADS-box TF which regulates the
growth of sepals, petals and floral meristems (Li et al. 2016). AP1/FUL
genes (also known as FUL-like genes)
in monocot can be divided into three groups, including FUL1, FUL2 and FUL3 (Wu et al. 2017), which tend to have
unnecessary functions in induction and differentiation of flower meristem. The
expression patterns of FUL-like genes
in monocot are more diverse than those of dicots, revealing that they may exert
different functions in monocot (Wu et al. 2017). We employed the
high-throughput sequencing to obtain a set of lncRNA
whose targeted gene was FUL1,
down-regulated in drought. Although most members of MADS-box TF show a positive
expression to several abiotic stresses like cold, salt and drought, some of them
present a negative regulatory relationship between response and stresses, such
as hormones. The expression of 12 HbMADS-box genes was up-regulated after drought stress in
rubber tree (Wei et al. 2018), while MADS-box gene ZMM7-L in maize down-regulated by exogenous ABA (Zhang et al. 2012).
Therefore, we speculated that there are some MADS-box genes playing negative or
reversed role in stress response. The WUSCHEL-related
homeobox
(WOX) TFs, which are specific in plant, play an essential part in a
variety of physiological and developmental processes (Li et al. 2019). It may
not only integrate auxin and cytokinin
signaling to regulate crown root development, but also respond to various abiotic
stresses, like drought, cold, and high salinity (Cheng et al. 2014). The expression of WOX11in rice could be quickly induced by
drought stress, and the up-regulated expression further enhance the drought
resistance, nevertheless the wox11
mutant is more sensitive to drought (Cheng et al. 2014). Cheng et al. (2016)
showed that, compared to wild type, the transgenic plants with overexpressed OsWOX11 significantly increased root hair length and amount, which
are conducive to water absorption and prompting drought resistance. Previous
research also indicated that, as a high hierarchical regulator, WOX11 not only involves in root
development and carbon metabolism, but also regulates other transcription
factors as well as genes in different pathways (Jiang et al. 2017). From the current
research, we screened a targeted mRNA named WOX11,
which responded the drought stress and expressed higher in B. distachyon.
According to the conservation of WOX genes, it might share the similar
function like in rice and A. thaliana.
However, the further research is required to verify its role in drought
resistance.
Ribosome proteins function in two ways in organisms: they not only
participate in protein synthesis as an important component of ribosomes; but
also take part in the growth and development processes (Lee et al. 2018). The
ribosomal protein genes highly expressed in rice plants infected with leaf
blight, which showed a high response to stress signals, indicating that
ribosomal proteins were a valuable resource to regulate plant stress resistance
(Moin et
al. 2016). Kang et al. (2016) isolated ribosomal P3 protein from Arabidopsis and proved that it can
protect cells from high and low temperature stress. Although the expression of
ribosomal proteins is relatively stable, environmental stress still affects its
transcription. In this study, the genes encoding ribosomal proteins S2, S18,
S19/S15, L33 and L23/25 are target genes for drought-responsive lncRNAs, but drought inhibited the expression of these
target genes, which indicates that these ribosomal proteins may play a negative
regulatory role under drought stress. However, the specific mechanisms and
approaches still remain to be elucidated.
Apart from photosynthetic systems, riboproteins,
MADS-box TFs and WOX TFs responding to drought stress in this experiment, some
enzymes involved in nitrogen metabolism and carbon metabolism also showed
distinct expression patterns. Caseinolytic protease (ClpP) is a serine protease with high conservation. The
expression diversity of ClpP
is particularly evident in higher plants, especially under stress conditions.
Overexpression of rice ClpD1 protein in transgenic Arabidopsis plants could increase their tolerance to salt and
drought stress (Mishra and Grover 2016).
However, under long-term drought stress in the field, the expression level of ClpP protein remained stable in drought-tolerant wheat
varieties Zlatitza and Yantar,
but decreased in drought-sensitive wheat cultivars Miziya
and Dobrudjanka (Vassileva et al. 2012). The target gene
encoding ClpP was down-regulated after drought
treatment in our experiment, similar to the results of Vassileva et al. (2012).
As a sessile organism, the complexity of the Clp
family in plant may be reflected in the response to multiple environmental
stresses during its long-term evolution (Williams et al. 2019). The glycoside
hydrolases in plant not only take part in the metabolism of various
carbohydrates, such as cell wall polysaccharide metabolism, but also involve in
the biosynthesis and regulation of glycans, adversity
defense, signal transduction, secondary metabolism and so on (Yue et al. 2019).
A glycoside hydrolase that responds to drought stress obtained in this
experiment is described as an alpha amylase (AMY1), which is down-regulated
after drought treatment but is inconsistent with previous studies. In terms of
previous reports, three genes (aAmy1,
aAmy3, and aAmy8) that encode α-amylase were up-regulated significantly
in rice under submergence
conditions (Du et al. 2014). It is speculated that the high expression of alpha amylase under
stress conditions may be related to stress response, which can provide energy
for plant growth by hydrolyzing starch into metabolizable
sugars to overcome subsequent stresses (Shaw et al. 2017). However, the
results of this study do not correspond to this speculation, suggesting that
alpha-amylase may have additional regulatory mechanisms under adverse
conditions.
Transcriptional regulation and post-transcriptional regulation are
important ways to regulate gene expression (Gonçalves et al. 2017). In recent years,
increasing genomics and genetic studies have focused on lncRNAs,
which are the essential part of transcriptional and post-transcriptional
regulation rather than transcriptional noise in plant. LncRNAs
not only participate in the growth and development of plant, but also involve
in the response of stresses (Di et al. 2014; Yuan et al. 2016). LncRNAs may be
involved in the regulation of plants at transcriptional and post-transcriptional
level through a variety of pathways. Some lncRNAs
possess a binding site for miRNAs, so they can act as
targets for specific miRNAs via target mimics, in order to prevent interactions between miRNAs and their targets (Wang et al. 2017a). For example, a lncRNA in rice, osa-eTM160, has been demonstrated as an
endogenous target mimics of osa-miR160 and represses expression of osa-miR160
to up regulate the osa-ARF18,
which is related with the reproductive development of rice (Wang et al. 2017b). A few plant lncRNAs could serve as precursors to small ncRNAs (sRNAs), such as miRNAs and siRNAs, to produce sRNA (Wang et al. 2017a).
Certain antisense lncRNAs in plant would respond to environmental stress by
interacting with sense mRNAs (Wang et al. 2017a).
Antisense RNA might take part in regulation of gene expression at various
levels, including genomic imprinting, transcription, transcriptional
interference, RNA stability, alternative splicing, and chromatin modification (Wight and Werner 2013; Zhang et al. 2013).
Wang et
al. (2014a) confirmed antisense lncRNAs in
A. thaliana in response to light,
which may form a complex with sense mRNA affecting gene expression in the
opposite strand. We also obtained some lncRNAs that
responded to drought in the leaves of B. distachyon under drought treatment, and their
expression patterns were quite different. Even if lncRNAs
predicted to be the same target gene, there may be an opposite expression
pattern. For example, the target gene described as WOX TF, induced by drought,
is regulated by 27 lncRNAs, 12 of which are
down-regulated and 15 are up-regulated (Table S7). Although these lncRNAs may work together on the same target gene, the
results indicate that lncRNAs might regulate target
gene expression through multiple pathways and participate in diverse biological
processes. The fine regulation mechanism and action mode still need further
research to explore.
Conclusion
A total of 90 lncRNAs with differential expression were screened in this
research. According to the functional annotations, the predicted target genes
were divided into five categories, including three photosynthesis-related
genes, two transcription factors, one peptide, five ribosomal proteins, and one
glycoside hydrolase. Most of the target genes were down-regulated under drought
treatment except the transcription factor wox11. All of the valuable
information provides a theoretical basis for further characterizing candidate
genes involved in drought stress. However, further studies of annotated and
novel transcripts are required to elucidate their function in adaptation
mechanisms under drought stress. In conclusion, our research improves new
perspectives into the defense mechanisms of drought response in B. distachyon,
and also improves the comprehensive understanding concerning the drought
tolerance as well as resistance in other temperate cereal.
Acknowledgements
This research was
financially supported by the National Key R&D Program
of China (Grant No. 2017YFD0301100), National
Natural Science Foundation of China (Grant No. 31401323), the Natural Science
Foundation of Henan Province (Grant No.
182300410040), the Key Research Projects of Higher Education in Henan Province
(Grant No. 18B210002), the Fund of Henan University of Science and Technology
(Grant Nos. 13480072 and 09001814), HAUST discipline improvement and promotion
plan A (Grant No. 13660002), and Twelfth Five-Year National Science and
Technology Pillar Programme (Grant No. 2011BAD16B07).
References
An
N, S Fan, Y Wang, L Zhang, C Gao, D Zhang, MingyuHan (2018). Genome-wide
identification, characterization and expression analysis of long non-coding RNAs in different
tissues of apple. Gene 666:44‒57
Blümel M, N Dally, C Jung (2015). Flowering time
regulation in crops – what did we learn from Arabidopsis? Curr Opin Biotechnol 32:121‒129
Cheng
S, DX Zhou, Y Zhao (2016). WUSCHEL-related
homeobox gene WOX11 increases rice
drought resistance by controlling root hair formation and root system
development. Plant Signal Behav 11;
Article e1130198
Florea L, L Song, SL Salzberg (2013).
Thousands of exon skipping events differentiate among splicing patterns in
sixteen human tissues. F1000Res 2;
Article 188
Gonçalves E, ZR Nakic, M Zampieri, O
Wagih, D Ochoa, U Sauer, P Beltrao, J Saez-Rodriguez (2017). Systematic
analysis of transcriptional and post-transcriptional regulation of metabolism
in yeast. PLoS Comput Biol 13; Article
e1005297
Lee CH, M Kiparaki, J Blanco, V Folgado,
Z Ji, A Kumar, G Rimesso, NE Baker (2018). A regulatory response to ribosomal
protein mutations controls translation, growth, and cell
competition. Dev Cell 46:456‒469
Matiu M, DP Ankerst, A Menzel (2017). Interactions
between temperature and drought in global and regional crop yield variability
during 1961‒2014. PLoS One 12;
Article e0178339
Sasaki T, BA Antonio (2004). Rice
Genome as A Model System for Cereals. Springer,
Dordrecht, The Netherlands
Scholthof KBG, S
Irigoyen, P Catalan, KK Mandadi (2018). Brachypodium: a monocot grass model genus for plant
biology. Plant Cell 30:1673‒1694
Wang
YQ, XD Fan, F Lin, GM He, W Terzaghi, DM Zhu, XW Deng (2014b). Arabidopsis noncoding RNA mediates
control of photomorphogenesis by red light. Proc
Natl Acad Sci USA 111:10359‒10364
Wang
M, W Zhao, L Gao, L Zhao (2018). Genome-wide profiling of long non-coding RNAs
from tomato and a comparison with mRNAs associated with the regulation of fruit
ripening. BMC Plant Biol 18; Article
75
Wang JJ, XW Meng, OB Dobrovolskaya, YL
Orlov, M Chen (2017a). Non-coding RNAs and their roles in stress response in
plants. Genomics Proteomics
Bioinform 15:301‒312
Wang M, HJ Wu, J
Fang, CC Chu, XJ Wang (2017b). A
long noncoding RNA involved in rice reproductive development by negatively
regulating osa-miR160. Sci Bull
62:470‒475
Wei B, RZ Zhang, JJ Guo, DM Liu, AL Li,
RC Fan, L Mao, XQ Zhang (2014). Genome-wide analysis of the mads-box gene
family in Brachypodium distachyon. PLoS One, 9; e0084781
Wei M, Y Wang, R
Pan, W Li (2018). Genome-wide
identification and characterization of MADS-box family genes related to floral
organ development and stress resistance in Hevea
brasiliensis Müll. Arg. Forests
9:1-16
Wight M, A Werner (2013). The functions
of natural antisense transcripts. Essays
Biochem 54:91‒101
Zanetti G, A Allverti (2018). Chemistry and Biochemistry of Flavoenzymes.
CRC Press, Boca Raton, Florida, USA